%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% O-Ring Production Networks
%   Author: Demir, Fieler, Xu, Yang
%   Last updated: January 2021
%
% Purpose: Endogenous Targeting (step 2)
%   - Select best fit model for each nuc
%   - Load saved initial equilibrium
%   - Solve baseline counterfactual
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all; diary off; clc; clear; 

    % load estimates (exogenous targeting)
    load('../../output/estimates/fullmodel_estimates.mat', 'est');            
 
    % load endogenous targeting calibration summary
    load('../../output/temp/endogtarget_initeqm_grid.mat')     
   
    
    
%% Initialize

    endog.nuc = zeros(Nc,1);
    endog.nuv = zeros(Nc,1);
    endog.fval = zeros(Nc,1);
    endog.cf1_wage_exp = zeros(Nc,1);
    endog.cf1_wage_nonexp = zeros(Nc,1);
    endog.cf1_wage_all = zeros(Nc,1);
    endog.cf1_qlevel_all = zeros(Nc,1);
    endog.eqm_tauindex = cell(Nc,1);
    endog.cf1_tauindex = cell(Nc,1);
    
    
%% Loop (Counterfactual)
for nc = 1:Nc    
tic    
    
    %% Pick best fit calibration
    
        % best fit nuv index
        fval_matrix_adj = fval_matrix + 100*(converg.iter_outer==100*ones(Nc,Nv)) + 100*(converg.iter_inner==500*ones(Nc,Nv));
        [~, nv] = min(fval_matrix_adj(nc,:));
        
        %nc=5;nv=4; cycle of 15 iterations
        
        % nu_c and nu_v
        nu_c = nuc_matrix(nc, nv);
        nu_v = nuv_matrix(nc, nv);
        fprintf('\nCounterfactual(%.0f,%.0f): nu_c = %.4f, nu_v = %.4f\n', nc, nv, nu_c, nu_v)
        fprintf('Saved calibration: fval = %.4f, err_outer = %e, err_inner = %e, iter_outer = %.0f, iter_inner = %.0f\n', ...
           fval_matrix(nc,nv), converg.err_outer(nc,nv), converg.err_inner(nc,nv), converg.iter_outer(nc,nv), converg.iter_inner(nc,nv));    
        
        
    %% Initial equilibrium
    fprintf('\nInitial Equilibrium\n');
    
        % parameters
        define_param
        define_param_est_endogtarget     
        param.nu_c = nu_c; %calibrate nu_c
        param.nu_v = nu_v; %calibrate nu_v

        % data moment
        define_data_moment

        % constants
        define_omega_dist 
        define_const_endogtarget 

        % solve initial equilibrium
        eqm = solve_eqm_outer_endogtarget(param, const);    
        define_wage_schedule
        const.phi_v = eqm.phi_v; %add phi_v to const

        % calculate model moments and objfcn
        model = cal_model_moment(param, const, eqm); 
        bartik = cal_chg_bartik(param, const, eqm);     
        fval = cal_objfcn(data, model, bartik);    
    
        % calculate firm vars
        firm = cal_firm_vars(param, const, eqm);

    
    %% Baseline counterfactual 
    fprintf('\nBaseline Counterfactual\n');

        % 5 percent export demand shock
        const.D_F = param.shock*const.D_F;    

        % solve the equilibrium
        eqm_cf1 = solve_eqm_cf1_outer_endogtarget(param, const, eqm);
        const.phi_v = eqm_cf1.phi_v; %update phi_v to const

        % calculate firm vars
        firm_cf1 = cal_firm_vars(param, const, eqm_cf1); 
        
        % counterfactual changes
        ge_cf1 = cal_chg_ge(eqm, firm, eqm_cf1, firm_cf1);       

    
    %% Store results
    
        % display
        fprintf('\nparam.nu_c = %.2f, param.nu_v = %.4f, fval = %.4f\nInitial equilibrium: n(qmax) = %e, err_outer = %e, err_inner = %e, iter_outer = %.0f, iter_inner = %.0f\n', ...
            param.nu_c, param.nu_v, fval, eqm.nq(param.Q), eqm.err_outer, eqm.err_inner, eqm.iter_outer, eqm.iter_inner);    
        fprintf('Baseline counterfactual: n(qmax) = %e, err_outer = %e, err_inner = %e, iter_outer = %.0f, iter_inner = %.0f\n', ...
            eqm_cf1.nq(param.Q), eqm_cf1.err_outer, eqm_cf1.err_inner, eqm_cf1.iter_outer, eqm_cf1.iter_inner);      
    
        % store
        endog.nuc(nc) = param.nu_c;
        endog.nuv(nc) = param.nu_v;
        endog.fval(nc) = fval;
        endog.cf1_wage_exp(nc) = ge_cf1.overall_wage_exp;
        endog.cf1_wage_nonexp(nc) = ge_cf1.overall_wage_nonexp_ctn;
        endog.cf1_wage_all(nc) = ge_cf1.overall_wage_all;
        endog.cf1_qlevel_all(nc) = ge_cf1.overall_qlevel_all; 
        endog.eqm_tauindex{nc} = eqm.tau_index;
        endog.cf1_tauindex{nc} = eqm_cf1.tau_index;
        
        
toc
end



%% Output figures
diary off

    % create folder
    mkdir('../../output/figures')    

    % output
    out_endogtarget_figures


